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1. Introduction 

Recent years have seen the introduction of powerful new tools for studying degenerate 
quantum gases. For example, on the experimental side correlation measurements offer 
a new experimental probe of many-body effects [Il[2l[3l|4l|5l[6l[7l[8l[9l [lOl [H] • On the 
theoretical side, the challenges posed by the new experimental techniques are being met 
by quantum dynamical simulations of large numbers of interacting particles in realistic 
parameter regimes. These are becoming possible due to the advances in computational 
power and improvements in numerical algorithms (for recent examples, see [T2| [T3 | [H]). 

In this paper we study metastable helium (^He*), which is currently unique in 
quantum atom optics in that it permits a comparison of experimentally measured [15] 
and theoretically calculated quantum correlations. This is one of the first examples in 
which experimental measurements can be considered in the context of first-principles 
calculations. Our goal in this paper is to confront a theoretical analysis with the results 
of recent experiments on atomic four-wave mixing via a collision of two Bose-Einstein 
condensates (BECs) of metastable ^He* atoms [15]. Figured] is a schematic momentum 
space diagram of these experiments. Two condensates, whose atoms have approximately 
equal but opposite momenta, ki and k2 ~ — ki, interact by four- wave mixing, while 
they spatially overlap, to produce correlated atomic pairs with approximately equal but 
opposite momenta, ka and k4, satisfying momentum conservation, ki -|- k2 = ks -|- k4. 
Figure [1] corresponds to the experimental data shown in figure 2 of [15], since after 
time-of-fiight expansion, atomic momentum is mapped into atomic position. 

We perform first-principles quantum simulations of the collision dynamics using 
the positive P-representation method [161 [IS [H]. The advantage of this method 
is that given the Hamiltonian of the interacting many-body system, no additional 
approximations are imposed to simulate the quantum dynamics governed by the 
Hamiltonian. The drawback on the other hand, is that it usually suffers from large 
sampling errors and the boundary term problem [19] as the simulation time tgim 
increases, eventually leading to diverging results. 

An empirically estimated upper bound for the positive-P simulation time (with 
controllable sampling error) of condensates with s-wave scattering interactions is given 
approximately by [20] 

tsim < 2.5m{AVY/y[A7rhapl% (1) 

where m is the atom mass, a is the s-wave scattering length, po is the condensate peak 
density, and AV = AxAyAz is the volume of the elementary cell of the computational 
lattice, with lattice spacings of Ax, Ay, and Az. Applying this formula to metastable 
helium, we see that this is a particularly challenging case among commonly condensed 
atoms due to its small atomic mass and relatively large scattering length. Our 
simulations are restricted to short interaction times (typically < 25 fis), which are 
about 6 times shorter than the experimental interaction time of p^. Despite this, 
our positive-P simulations provide useful insights into the experimental observations; 
in addition, they can serve as benchmarks for approximate theoretical methods (such 
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Figure 1. Schematic momentum space diagram of the atomic four- wave mixing 
interaction. Optical Raman pulses generate untrapped condensates with momenta 
ki and k2 — — ki parallel to the x-axis (dark disks). These undergo a four-wave 
mixing interaction to produce correlated atomic pairs on a spherical shell of radius ki . 

as the Hartree-Fock-Bogoliubov method [22l [23], [21]) to estabhsh the range of their 
vahdity. 

We calculate atom pair correlations within the scattering halo produced 
spontaneously during the collision (see figure [I]). The scattering halo is a spherical 
shell in momentum space. In the limit of small occupation of the scattered modes, the 
s-wave nature of the collisions ensures an approximately uniform atom population over 
the halo. We consider the strength and the width of the correlation signal, as well as 
the momentum width of the halo. We also analyze relative atom number squeezing and 
the violation of the classical Cauchy-Schwartz inequality. 

In Sec. 2 of this paper we will summarize the experimental results we wish to 
analyze. In Sec. 3 we discuss order of magnitude estimates. In Sec. 4 we describe 
simulations using the positive P-representation method, and in Sec. 5 we discuss the 
results of our simulations. Sec. 6 summarizes our work. 

2. Summary of experimental results 
2.1. Overview of the experiment 

The starting point of the experiment is a '^He* condensate of 10^ to 10^ atoms confined in 
a magnetic trap whose frequencies are: uJx/27t = 47 Hz and ujy/2'K = uJz/2tt = 1150 Hz. 
A sudden Raman outcoupling drives the trapped ^He* from the = 1 Zeeman sublevel 
into the magnetic field insensitive state = 0. [I5]. The Raman transition also splits 
the initial {m^ = 1) condensate into two roughly equally populated condensates with 
opposite velocities along the x direction. The magnitude of each velocity is equal to the 
recoil velocity Vr = 9.2 cm/s, defined by the momentum of the photons used to create 
the colliding condensates hkr, kj. = 5.8 x 10^ m~^. The relative velocity 2vr of the two 
condensates is about 8 times higher than the speed of sound Cs = of the initial 

condensate, ensuring that elementary excitations of the condensates correspond to free 
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particles. 

During the separation of the condensates, elastic collisions occurring between atoms 
with opposite velocities scatter a small fraction (5%) of the total initial atom number 
into the halo. The system is shown in three-dimensions in an accompanying video of the 
experimental results after a 320 ms time of flight 0. For the purposes of this paper, the 
experiment consists in the acquisition of the three dimensional positions of the particles 
scattered into the collision halo after the time of flight. This information permits the 
reconstruction of the 3D momentum vectors of the individual particles after they have 
ceased interacting with each other. 



2.2. Main results 

Knowledge of the momentum vectors in turn permits the construction of two-particle 
correlation functions in momentum space. The correlation function shows features for 
particles traveling both back to back and collinearly. The back-to-back correlation 
results from binary, elastic collisions between atoms, whereas the collinear correlation 
is a two particle interference effect involving members of two different pairs: a Hanbury 
Brown-Twiss correlation [25]. Both correlation functions are anisotropic because of the 
anisotropy of the initial colliding condensates. 

To quantify these correlations, we first introduce the unnormalized normally- 
ordered second-order correlation function between the densities at two points in 
momentum space, 

G'(2)(ki,k2) = (:n(ki)n(k2):). (2) 

Here, n(k) = a^^(k)a(k) is the momentum density operator, d^(k) are a(k) are the 
Fourier transforms of the field creation and annihilation operators ^^^(x) and ^(x), and 
the colons :: stand for normal ordering of the operators according to which all creation 
operators stand to the left of the annihilation operators, so that 

(: n(ki)n(k2) :) = (at(ki)at(k2)a(k2)a(ki)). (3) 

Because of a low data rate, the correlation measurements must be averaged over the 
entire collision sphere to get statistically significant results. The average collinear (CL) 
second-order correlation as a function of the relative displacement Afcj in the fcj-direction 
{i = X, y, z) is defined as 

Jd^k G'(2)(k,k + eiAfc,) 
= /^3k (n(k))(n(k + e,AA;,))' 

V 

I A 3 dimensional, animated rendition of the atomic positions 320 ms after release from the trap. The 
vertical positions are derived from the arrival times as described in [15j . Each point corresponds to 
the detection of one atom and the animation shows the sum of 50 separate runs. The ellipsoids at the 
sides are the colliding condensates. The ellipsoids at the top and bottom result from imperfect Raman 
polarizations and stimulated atomic 4 wave mixing (see [15j). The 4 condensates are excluded from 
the analysis given in the text. 
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Figure 2. Cross section of tlie scattering halo. A sloped background is present due 
to thermal atoms in the trap. This background has been fit to a straight line and 
subtracted in order to estimate the rms width, Sk ~ 0.067A:,.. 



(2) 

where e, is the unit vector in the ki direction. The normahzation of g^j^i^Aki) ensures 
that for uncorrelated densities g(~;[{Aki) = 1. The integration domain D in (jl]) selects 
a certain region of interest in momentum space and can be defined, for example, to 
contain only a narrow spherical shell and to eliminate the initial colliding condensates. 
Due to the averaging, the dependence of the correlation functions on the direction k is 
lost. 

(2) 

The average back-to-back (BB) correlation function g^^{Aki) between two 
diametrically opposite points, one of which is additionally displaced by Aki in the ki- 

(2) 

direction, is defined similarly to g^^^Aki): 

Jd^k G^^\k, -k + eiAki) 



/#k {h{k)){h{ 

V 



-k + BiAki))' 



The experimental observations can be summarized as follows. The width of both 
correlation functions along the axial direction of the condensate, the x-axis, is limited by 
the resolution of the detector and hence contains little information about the collision. 
In the radial direction (with respect to the symmetry of the colliding condensates), 
one observes a peak which can be fitted to a Gaussian function with rms widths u^f' 
and (T^f for collinear and back-to-back cases respectively. The experimental results are 
summarized in the following table 





<'z/kr 


^CL, BB 

^ y,z / y,z 


0.081 ±0.004 


0.069 ±0.008 


0.85 ±0.15 



(6) 



One can also use the data to deduce the averaged radial width 6k of the scattering 
halo. Figure [2] shows a cross section of the halo, averaged over all accessible scattering 
angles. The presence of the unscattered condensates prevents observation of the shell 
along the x-axis, but along the accessible directions we find 6k 



0.067/c^. 
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3. Qualitative analysis 

In this section we discuss some simple, analytical estimates of the measured quantities. 
In later sections we will do more precise, numerical calculations which will verify the 
results of this section. 



3.1. Width of the pair correlation functions 



As discussed in [15j, the width of the back-to-back and coUinear correlation functions 
should be on the order of the momentum width of the initial condensate, which 
in turn is proportional to the inverse width of its spatial profile. For a Gaussian 
density profile of the initial condensate in position space p(x) = (\l/'''(x)\E'(x)) = 
Po exp[— ^j^j, y 2 rf/{2wf)], and therefore a Gaussian density distribution in momentum 
space, n{k) = {n(k)) oc Gw[^J2i=xyz^i/('^^i)]' with (Xj = l/wi, an approximate 
theoretical treatment based on a simple ansatz for the pair wavefunction predicts a 
Gaussian dependence of the back-to-back (BB) and collinear (CL) correlation functions 
on the relative wavevectors Ah [25]: 



^^^^(k, -k + UiAki) oc exp 
G^^\k,k + niAki 



Akl 



oc exp 



2(af^)2 

CL\2 



2 a. 



(7) 
(8) 



The widths of the back-to-back {cf^) and collinear (crf^) correlations are related 
to the momentum-space width ai of the initial (source) condensate via [25] 

a^/a, = V2, (9) 



a. 



CL 



(10) 



and therefore the width of the back-to-back correlation is times smaller than the 
width of the collinear correlation. 

In Sec. 15.11 the initial momentum-space widths are found to be = 0.0025/Cr 



and 0", 



y,z 



0.055/cr, assuming 



measured widths in units of cr^^^, we can rewrite 



9.84 X 10^ atoms. Expressing the experimentally 

as 





<^iz/(^y,z 


^CL, BB 

y,z / y,z 


1.47 ±0.07 


1.25 ±0.15 


0.85 ±0.15 



:iii 



and therefore, ([9]) is in agreement with the measured width of the radial back-to-back 
correlation function, whereas flTUl) overestimates the width of the collinear correlation 
function by almost 60%. As we show below, first-principles simulations using the 
positive-P method and incorporating atom-atom interactions result in widths which 
are closer to the experimental values. 

The discrepancy between the two theoretical approaches (which apparently is larger 
for the collinear correlations than for the back-to-back ones) comes mostly from the fact 
that the above calculation is made for a Gaussian shape of the initial BEG density 
profile, whereas in practice and in the positive-P simulations the spatial density of a 
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harmonically trapped condensate is closer to an inverted parabola (as in the Thomas- 
Fermi limit) rather than to a Gaussian. An alternative theoretical model [26], based on 
the undepleted source condensate approximation and a numerical solution to the linear 
operator equations of motion for scattered atoms, also confirms that for short times 
the momentum-space correlation widths are narrower if the source condensate has a 
parabolic spatial density profile, compared to the case of a Gaussian density profile. 



3.2. Width of the scattered halo 

A second, experimentally accessible quantity in a BEG collision is the width hSki in 
momentum space of the halo on which the scattered atoms are found. Glearly the 
momentum spread ai {in i = x, y or z direction) of the colliding condensates imposes a 
minimum width 

6h > cTi. (12) 

This hmit suggests that the halo could be anisotropic. As noted above however, the 
experiment in ^15] is not highly sensitive to such an anisotropy, and measures the width 
chiefly in the y, ^-directions. 

Other physical considerations also affect this width, and suggest that the halo 
should rather be isotropic, in which case we can drop the index from 6k. Here we 
discuss two mechanisms that impose a finite radial width on the halo. 

If the pairs are produced during a finite time interval At, the total energy of the 
pair is necessarily broadened by h/At. This is true even if the relative momentum is 
well defined. For a mean fc-vector kr, the finite interaction time between the colliding 
BEGs results in a broadening of 

rfi 

hkrAt ' ^ ^ 

where we assumed 5k /kr <^ 1. In the experiment, the collision time is sufficiently 
long that the above effect does not impose a limitation on the width of the sphere. In 
the positive-P simulations however, numerical stability problems limit the maximum 
collision time that can be simulated, as discussed in Sec. [5], and this time does indeed 
impose a width on the halo. For short collision times, where the scattering is in the 
spontaneous regime, our numerical results for the width 5k are in good agreement with 
the simple estimate of Equation f|T3|) . 

For long collision times it can happen that so many atoms are scattered that Bose 
enhancement and stimulated effects become important. In this case the width of the 
scattering shell can be estimated by a slightly more involved approximate approach 
based on analytic solutions for the uniform system within the undepleted "pump" 
(source condensate) approximation [27j. Under this approximation, the present system 
is equivalent to the dissociation of a condensate of molecular dimers studied in detail 
in [T3l [28l [29] . The latter system in turn is analogous to parametric down-conversion in 
optics [Sn]. The details of the approximate solutions, common to condensate collisions 
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and molecular dissociation, and the relationship between them are given in Appendix C 
The resulting width of the halo found from this approach is 

6k ^ (14) 

We see that in this regime, the width is proportional to the scattering length a and the 
peak density po, but it no longer depends on the collision duration. 

The physical interpretation of Equation (IT^ is that with the stronger effective 
coupling (or nonlinearity) apo, one can excite and amplify spectral components that 
are further detuned from the exact resonance condition hA^ = (or further "phase 
mismatched"). The inverse dependence on collision momentum kr can be understood 
via the quadratic dependence of the energy on momentum k: to get the same excitation 
at a given energy offset hA^, (10.30 . one requires smaller absolute momentum offset 6k 
at larger kr than at small kr. 

Positive-P simulations covering the transition from the spontaneous to stimulated 
regimes are available for ^^Na condensate collisions as in [T3]. The numerical results 
in this case are in agreement with the simple analytic estimate of Equation iHM . 
More specifically, we find that for collision durations between 300 /is and 640 fis the 
actual numerical results for the width of the spherical halo vary, respectively, between 
6k/kr ^ 0.13 and 6k/kr ~ 0.087, whereas Equation f|T4l) predicts 6k/kr ~ 0.096. 

For ^He*, on the other hand, the small mass and the larger scattering length of 
^He* atoms limit the maximum simulation time to tsim ^ 25 /is. This is far from the 
stimulated regime and therefore we do not have a direct comparison of the numerical 
results with Equation (IHj). The experiment is also not in the stimulated regime. We are 
nevertheless tempted by the numerical ^^Na result to extrapolate Equation f|T^ to ^He* 
BEG collisions in the long time limit and we obtain 6k/ kr — 0.05. Adding this width 
in quadrature to the momentum width of the initial condensate, ay^z — 0.055kr, gives 
^^/(0^05k/)^~+J0^055kr^ = 0.074^^, not far from the experimentally observed radial 
momentum width of 6k ^ 0.067/cr- We thus suggest that the mechanism leading to 
Equation f|T^ may play a role in the experiment. 

4. Model 

The effective field theory Hamiltonian governing the dynamics of the collision of BECs 
is given by 

H = J dx jl^lV^p + :^^t^t^^| ^ (15) 

where ^(x, t) is the field operator with the usual commutation relation 
[^(x, t), ^^(x', t)] = ^'•^''(x — x'), m is the atomic mass, the first term is the kinetic 
energy, and the second term describes the s-wave scattering interactions between the 
atoms. The trapping potential for preparing the initial condensate before the collision is 
omitted since we are only modeling the dynamics of the outcoupled condensates in free 
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space. The use of the effective delta function interaction potential U{x — y) = UoS{x — y) 
assumes a UV momentum cutoff fc™**^. In our numerical simulations the momentum cut- 
off is imposed explicitly via the finite computational lattice. If the lattice spacings (Ax, 
Ay, Az) in each spatial dimensions are chosen to be much larger than the s-wave scat- 
tering length a, then the respective momentum cutoffs satisfy k^y% <^ 1/a. In this case 
the coupling constant Uq is given by the familiar expression Uq ~ Anha/m [31] without 
the need for explicit renormalization. 

To model the dynamics of quantum fields describing the collision of two BECs, 
we use the positive P-representation approach [16]. In this approach the quantum field 
operators ^'(x, t) and ^^(x, t) are represented by two complex stochastic c-number fields 
\E'(x, t) and \E'(x, t) whose dynamics is governed by the following stochastic differential 
equations [14j : 

= ^V^vl/ - ^U,^J^J^> + ^fllW^ Ci(x, t), (16a) 

^^^= -7^V2^ + ^f/oW + yS^^C2(x,t). (166) 
ot 2m ^ 

Here, Ci(x, t) and C2(x, t) are real independent noise sources with zero mean, (Cj(x, t)) = 

0, and the following nonzero correlation: 

(0(x,t)a(x',t')) = - x')5(t - t'). (17) 

The stochastic fields \E'(x, t) and \E'(x, t) are independent of each other [\I'(x, t) ^ 
^*(x, t)] except in the mean, (\E'(x, t)) = (\E'*(x, t)), where the brackets (...) refer to 
stochastic averages with respect to the positive P-distribution function. In numerical 
realizations, this is represented by an ensemble average over a large number of stochastic 
realizations (trajectories). Observables described by quantum mechanical ensemble 
averages over normally-ordered operator products have an exact correspondence with 
stochastic averages over the fields \E'(x, t) and \E'(x, t): 

([^t(x,t)]-[^(x',t)]") = ([^(x,t)]-[vl/(x',t)]"). (18) 

The initial condition for our simulations is a coherent state of a trapped condensate, 
modulated with a standing wave that imparts initial momenta ±kr (where kr = mVr/h 
and Vr is the collision velocity) in the x direction, 

M/(x, 0) = (^(x, 0)) = ^/M^ (e*'^^ + e"^'^^) , (19) 

with \^(x, 0) = \E'*(x, 0). Here, po(x) is the density profile given by the ground 
state solution to the Gross-Pitaevskii equation in imaginary time. The above initial 
condition models a sudden Raman outcoupling of a BEG of trapped ^He* atoms in the 
= 1 sublevel into the magnetic field insensitive state m^, = 0, using two horizontally 
counter-propagating lasers and a third vertical laser [15]. In this geometry, the Raman 
transitions split the initial {rrix = 1) condensate into two equally populated condensates 
and simultaneously impart velocities of ±Vr onto the two halves. As a result the two 
outcoupled condensates undergo a collision and expand in free space. Accordingly, in 
our dynamical simulations, the field \l/(x, t) represents the atoms in the untrapped state 
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rrirc = 0, having the s-wave scattering length of aoo = 5.3 nm ([15] and references 
therein), while the initial density profile po(x) refers to that of the trapped atoms 
in the rrix = 1 state having the scattering length of an = 7.51 nm [32]. The same 
distinction in terms of the scattering length in question applies to the definition of the 
interaction strength Uq ~ ATcha/m, in which a has to be understood as an for the 
trapped condensate or as aoo for the outcoupled cloud. 

In our simulations we assume for simplicity that the outcoupling from the trapped 
rrix = 1 state is 100% efficient, in which case the entire population is transferred into 
the = state and therefore we have to only model s-wave scattering interactions 
between the atoms in the rrix = state. In the experiment, on the other hand, the 
transfer efficiency is only about 60% and therefore the collisions between the atoms 
in the m^. = and m^. = 1 are not completely negligible and maybe responsible for 
some of the deviations between the present theoretical results and the experimental 
observations. 



5. Results and discussion 

5.1. Main numerical example 

Here we present the results of positive-P numerical simulations of collisions of two 
condensates of ^He* atoms (m ~ 6.65 x 10^^^ kg) as in the experiment of [I5]. The 
key parameters in our main numerical example are the collision velocity, Vr = 9.2 cm/s, 
and the peak density of the initial trapped condensate, po = 2.5 x 10^^ m^'^. The 
trap frequencies are matched exactly with the experimental values, u!x/2tt = 47 Hz and 
ujy /2n = uJz/271 = 1150 Hz. The s-wave scattering length for the magnetically trapped 
atoms in the = 1 sublevel is an = 7.5 nm; the s-wave scattering length for the 
outcoupled atoms in the = sublevel is aoo = 5.3 nm. Other simulation parameters 



are given in Appendix Appendix D 



The initial state of the trapped condensate is found via the solution of the Gross- 
Pitaevskii equation in imaginary time. Given the above trap frequencies and the peak 
density as a target, we find that the total number of atoms in the main example is 
N = 9.84 X 10'^. With these parameters, the average kinetic energy of colliding atoms is 
Ekin/ks = mv1/2kB — 2.0 x 10~^ K, which is about 7.4 times larger than the mean-field 
energy of the initial condensate EMp/ks = ^T^h^anpQ/mkB — 2.7 x 10"'' K. 

The duration of simulation in the main example is = 25 ps. This is considerably 
smaller than the estimated duration of collision in the experiment, 140 /xs (see 



Appendix A). The number of scattered atoms in our numerically simulated example 
at tf = 25 /is is ~ 1750, representing ~ 1.8% of the total number of atoms in the initial 
BEG. Operationally, the fraction of scattered atoms is determined as the total number 
of atoms contained within the scattering halo (see figure [3] showing two orthogonal slices 
through the momentum density distribution) after eliminating the regions of momentum 
space occupied by the two colliding condensates. We implement the elimination by 
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Figure 3. Slices through — (a) and fc^: = (b) of the 3D atomic density 
distribution in momentum space n(k,tf) after tf = 25 /iS coUision time. Due to the 
symmetry in the transverse direction (orthogonal to x), the average density through 
ky = coincides with that of kz = 0. The color scale is chosen to clearly show the 
halo of spontaneously scattered atoms and cuts off the high-density peaks of the two 
colliding condensates (shown in white on the left panel). 




0.5 1 1.5 

k/k 

r 



Figure 4. Angle averaged (radial) momentum distribution n{k) of the scattered atoms 
(solid line) and a simple Gaussian fit (dashed line) used to define the radial width 
6k = O.lOfcr of the halo around the peak momentum fco = 0.98fcr (see text). 



simply discarding the data points corresponding to {k^l > 0.99kr, which fully contain 
the colliding condensates. This cuts off a small fraction of the scattered atoms as well, 
but the procedure is simple to implement operationally and is unambiguous. 

In order to compare our calculated fraction of scattered atoms at t / = 25 /is with the 
experimentally measured fraction of 5% at the end of collision at ~ 140 fis, we first note 
that these time scales are relatively short and correspond to the regime of spontaneous 
scattering. The number of scattered atoms increases approximately linearly with time, 
therefore our calculated fraction of 1.8% can be extrapolated to about 10% to correspond 
to the expected fraction at ~ 140 fis. Next, one has to scale this value by a factor 0.6^ 
to account for the fact that in the experiment only 60% of the initial atom number 
was transferred to the rrix = state of the colliding condensates. Accordingly, our 
theoretical estimate of 10% should be proportionally scaled down to 4% conversion, in 
good agreement with the experimentally estimated fraction of 5% (see also [Appendix A]). 
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Ak/k Ak /k 

X r y,z r 

Figure 5. Back-to-back (BB) and collinear (CL) atom-atom pair correlation, 
dBB/CL^^^i) — 1 as a function of tlie displacement Aki {i = x,y,z) in units of 
the collision momentum fc^, after tf = 25 iJS collision time. The circles are the 
numerical results, angle-averaged over the halo of scattered atoms after elimination 
of the regions occupied by the two colliding condensates; the solid lines are simple 
Gaussian fits to guide the eye (see text). For comparison, we also plot the initial 
momentum distribution no(fci) of the colliding condensates; the actual data points are 
shown by the squares and are fitted by a dashed-line Gaussian. 



In figure m we plot the radial momentum distribution of scattered atoms (solid line), 
obtained after angle averaging of the full 3D distribution within the region < O.Skr- 
The numerical result is fitted with a Gaussian oc exp[— (A; — koY /{2Sk'^)] (dashed line), 
centered at ko = 0.98kr and having the radial width of 6k = O.IOA;^ — 5.8 x 10^ m~^, 
where k = |k|. The fitted radial width of 6k = O.lOfcr of the numerical simulation 
is in reasonable agreement with the simple estimate of Equation f|T3l) . which gives 
6k ~ 0.075A;^ for At = 25 /is. 

Figure [5] shows the numerical results for the back-to-back and collinear correlations 
(solid lines with circles), defined in Equations (jlj) and ([5]). Due to the symmetry of the y 
and z directions, the results in these directions are practically the same. In order to verify 
the hypothesis that the shape and therefore the width of the pair correlation functions 
is governed by the width of the momentum distribution of the source condensate, we 
also plot the actual initial momentum distributions of the source condensate in the two 
orthogonal directions (with the understanding that the horizontal axis Aki now refers 
to the actual wave- vector component fcj). The actual data points for the correlation 
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functions and for the momentum distribution of the source are shown by the circles and 
squares, respectively, and are fitted with Gaussian curves for simplicity and to guide the 
eye. The Gaussian fits for the correlation functions (solid lines) give: 

g^^liAh) - 1 = 9.2exp{-Afcf/[2(aD2]}, (20) 
ggliAk,) - 1 = exp{-AA:,V[2(aP)^]}, (21) 

where the correlation widths af^ and af^ are shown in the table fl22|) below. The 
Gaussian fits (dashed lines) for the slices of the initial momentum distribution no(fcj) oc 
exp{— A;,?/[2(crj)^]} are scaled to the same peak value as g^B/CLi^^ ~ ^ have 
ax = 0.0025/cr and ay^z = 0.055^^- 

By comparing the solid and the dashed lines we see that the shape of the correlation 
functions indeed closely follow the shape of momentum distribution of the source. More 
specifically, we find that the following results provide the best fit to our numerical data: 











1.18 


1.39 


1.27 


1.57 



(22) 



The ratios between the collinear and back-to-back correlation widths are a'^^ ja^^ ~ 
1.08 and c^f'/c^f — 1.13. The errors due to stochastic sampling on all quoted values 
of the correlation widths are smaller than 3%. 

The values for CTy^/ay^z and (yy^joy^z can be compared with the respective 
experimentally measured values of table flTT]) and we see reasonably good agreement, 
even though the numerical data are for a much shorter collision time. The remaining 
discrepancy between the numerical data at = 25 /is and the experimentally measured 
values after a ~ 140 /is interaction time may be due to the evolution of the condensates 
past 25 /is, not attainable within the positive-P method. The above numerical results for 
the correlation widths can be also compared with the simple analytic estimate based on 
the Gaussian Ansatz treatment of Equations ([9]) and (fTOl) . We find that the approximate 
analytic results overestimate the back-to-back and collinear widths by ~ 20% and 40%, 
respectively, in the present example. 

The amplitude of the correlation functions can also be inferred by simple models. 
In fact, the collinear correlation function is a manifestation of the Hanbury Brown and 
Twiss effect since it involves pairs from two independent spontaneous scattering events 

/OA 

and we expect an amplitude of S'cmO) ~ 2 [25]. This is in agreement with the positive- 
P simulations. The back-to-back correlation amplitude, on the other hand, can be 
substantially higher and display super- bunching (^)jIj(O) > 1) [131 [22] since the origin 
of this correlation is a simultaneous creation of a pair of particles in a single scattering 
event. 

In a simple qualitative model [15], the amplitude of the back-to-back correlation 
can be linked to the inverse population of the atomic modes on the halo. As we show 



m 



Appendix B this model follows the trends observed in our first-principles numerical 



simulations. 
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Figure 6. Same as in figure [3l except for tf = 12.5 /iS collision time. 




Figure 7. Same as in figure IH except for tf = 12.5 /xs collision time. The width and 
the peak of the fitted Gaussian here are: 6k = 0.20fcr and ko = 0.95fcr. 



5.2. Shorter collision time 

Here we present the results of numerical simulation for the same parameters as in our 
main numerical example from Sec. 15.11 except that the data are analyzed at tf = 12.5 
/is, which is half the previous interaction time. We found in Sec. 15.11 that ay^,ay^ 
and the width of the halo 5k are all nearly the same. In Sec. 3, however, we argue 
that the widths of the correlation functions and the halo are governed by different 
limits [Equations f|T0|) .fl9|) and f|T3l) or f|T^ . respectively]. The example in this section 
illustrates this point. 

Figure [6] shows two orthogonal slices of the s-wave scattering sphere in momentum 
space (cf. figure [3]), whereas figure [7] is the corresponding radial distribution after angle 
averaging. The most obvious feature of the distribution is that it is broader than at 
tf = 25 /xs and the fitted Gaussian gives the radial width of Sk = 0.20/Cj,. This is 
precisely twice the width in Figure H] and is in agreement with the simple qualitative 
estimate of Equation (fT3|) . 

The back-to-back and collinear correlation functions after tf = 12.5 fis collision 
time are qualitatively very similar to those shown in figure O except that the Gaussian 
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fits are 

g^PsiAh) - 1 = 35.6exp{-AA;,V[2(af^)2]}, (23) 
ggliAh) - 1 = exp{-AA;,V[2(aP)^]}, (24) 

with the correlation widths given by 











1.16 


1.28 


1.27 


1.48 



(25) 



The ratios between the widths are cr^^ /a^^ ~ 1.09 and CTy z/o'y^f — 1-16. 

For the correlation functions, the main difference compared to the case for 25 /zs 
is that the peak value of the back-to-back correlation is now larger, reflecting the lower 
atomic density on the scattering halo. The correlation widths, on the other hand, are 
practically unchanged, at least within the numerical sampling errors of the positive- 
P simulations; the errors are at the level of the third significant digit in the quoted 
values, which we suppress. The number of scattered atoms in this example is about 
850, which is approximately half the number at 25 /is, confirming the approximately 
linear dependence on time in the spontaneous scattering regime. 

5.3. Smaller collision velocity 

In this example, we present the results of simulations in which the collision velocity 
is smaller by a factor than before, v'^, = 6.5 cm/s {k'^ = 4.1 x 10^ ni~^), while all 
other parameters are unchanged. In practice, this can be achieved by changing the 
propagation directions of the Raman lasers that outcouple the atoms from the trapped 
state. As in the previous example, the halo width illustrates Equation (fT3l) . 

The results of positive-P simulations for the momentum density distribution at 
tf = 25 fis are shown in figures [8] and [91 The most obvious feature of the distribution is 
again the fact that it is now broader than in our main example of Sec. 15. 1[ The width of 
the Gaussian function fitted to the numerically calculated radial momentum distribution 
is given by 6k = 0.21k[.. This is again in excellent agreement with the simple analytic 
estimate of Equation f|T3|) . which predicts the broadening to be inversely proportional 
to the collision velocity. We also note that the peak momentum (relative to k'^) in the 
present example is slightly shifted towards the centre of the halo, ko = 0.92/c^, which is a 
feature predicted in [2^ to occur when the ratio of the kinetic energy to the interaction 
energy per particle is reduced. 

The back-to-back and collinear correlation functions in this example are again 
qualitatively very similar to those shown in figure [5l except that the Gaussian fits are 

g^Slm) - 1 = 9.0exp{-Ak^/[2iarm, (26) 
^gi(Afc,) - 1 = exp{-Ak^/[2{af^f]}, (27) 

with the correlation widths given by 









<^y,z/(^y,z 


1.16 


1.35 


1.31 


1.51 



(28) 
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Figure 8. Same as in figure[3l except for \/2 times smaller collision velocity, v'^. — 6.46 
cm/s (fc,(, — 4.09 X 10^ ni^^). The axis for the momentum components ki (i = x,y,z) 
are in units of the smaller recoil momentum than in figure [3l and therefore the absolute 
radius of the s-wave scattering sphere is smaller in the present example. 




r 



Figure 9. Same as in figure S] except for \/2 times smaller collision velocity 
v'^ (fc^ = 4.1 X 10^ m~^). The width and the peak of the fitted Gaussian are 
5k ^ 0.21fc; = 8.6 X 10^ ni-i; fco = 0.92fc;. 



where a^/k'^ — 0.0035 and a^/k'^ — 0.078. The ratios between the colhnear and back- 
to-back correlation widths are o'^^ ja^^ ^ 1.13 and o^^joy^ ~ 1.12. 

As we see from these results, the absolute widths of the correlation functions are 
practically unchanged compared to the main numerical example ( !22|) . This provides 
further evidence that, at least for short collision times, the correlation widths are 
governed by the momentum width of the source condensate, which is unchanged in 
the present example compared to the case of Sec. 15.11 

The number of scattered atoms in this example is about 1270, which is 
approximately -\/2 times smaller than in Sec. 15. H and corresponds to ~ 1.3% conversion. 
This scaling is in agreement with the rate equation approach [22], according to which 
the number of scattered atoms is proportional to the square root of the collision energy 
and hence to the collision momentum, which is times smaller here. 
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Figure 10. Same as in figure [3] except for the scattering lengtlis of ooo — 2.65 nm and 
ail = 3.75 nm, which are twice smaller than before. 



5.J^-. Smaller scattering length 

Finally, we present the results of numerical simulations for the same parameters as in 
our main numerical example from Sec. 15.11 except that the scattering lengths an and 
Ooo are artificially halved, i.e. aoo = 2.65 nm and an = 3.75 nm. The trap frequencies 
are unchanged and we modify the chemical potential to arrive at the same peak density 
of the initial BEC in the trap, po — 2.5 x 10^^ m~^. The total number of atoms is now 
smaller, ~ 3.5 x 10'*. One effect of changing the scattering length is that it changes 
the size and shape of the trapped cloud, and therefore also its momentum distribution. 
The shape is slightly closer to a Gaussian and therefore also to the treatment in [25] . 

Due to the smaller scattering length, the density distribution in position space of the 
initial trapped condensate is now narrower and conversely the momentum distribution 
of the colliding condensates is broader. On the other hand, the width of the halo 
(see figures [10 and [ID at tf = 25 /is) of scattered atoms is practically unchanged 
compared to the example of figure [H as it is governed by the energy-time uncertainty 
consideration (|T3l) . for the spontaneous scattering regime. The only quantitative 
difference is the lower peak density on the scattering sphere, which is due to the weaker 
strength of atom- atom interactions resulting in a slower scattering rate. The number 
of scattered atoms at 25 fis is ~ 180, corresponding to 0.51% conversion of the initial 
total number ~ 3.5 x 10^. The fraction of 0.51% itself corresponds approximately to 
a scaling law of ~ a^^'^, which is the same as the scaling of the total initial number of 
trapped atoms in the Thomas-Fermi limit for a fixed peak density. 

Since the widths of the correlation functions are governed by the width of the 
momentum distribution of the initial colliding condensates, we expect corresponding 
broadening of the correlation functions as well (see figure [T2!) . To quantify this effect, 
we fit the momentum distribution of the initial BEC by a Gaussian oc exp{— /c?/[2(crj)^]}, 
where = 0.0036^^ and ay^z = 0.068^^ (cf. with = 0.0025^^ and ay^z = 0.055kr 
in figure [5l which are ~ smaller). The Gaussian fits to the correlation functions in 



Atomic four-wave mixing via condensate collisions 



18 




0.5 1 1.5 

k/k 

r 



Figure 11. Same as in figure |4] except for twice smaller values of the scattering 
lengthes aoo and an. The width and the peak of the fitted Gaussian are 6k — O.lOkr 
and fco = 0.98fcr, which are the same as in figure [4] 



figure [12] are 

g^^lm) - 1 = 49exp{AA:,V[2(af«)2]}, 
g^^liAk,) - 1 = 0.94exp{Afc,V[2(aP)^]}, 

where the widths aj^^ and erf ^ are given by 











1.18 


1.53 


1.42 


1.81 



(29) 
(30) 



(31) 



We see that the relative widths are practically unchanged, implying that the 
absolute widths are broadened. The ratios between the collinear and back-to-back 
correlation widths are slightly increased and are given by a^^/cxf^ ~ 1.20 and 
- 1-18. 

These numerical results make the present example - with the diminished role 
of atom-atom interactions - somewhat closer to the simple analytic predictions of 
Equations ([9]) and (fTOl) based on a Gaussian ansatz for noninteracting condensates. 



5.5. Relative number squeezing and violation of Cauchy-Schwartz inequality 

Another useful measure of atom-atom correlations is the normalized variance of the 
relative number fiuctuations between atom numbers Ni and Nj in a pair of counting 
volume elements denoted via i and j, 



([A(iV,-iV,)]2) _^ , (:[A(iV,-iV,)p:) 



(32) 



{m + {N,) (AT,) + (iV,) 

where AX = X — {X) is the fluctuation. This deflnition uses the conventional 
normalization with respect to the shot-noise level characteristic of Poissonian statistics, 
such as for a coherent state, (Ni) + (Ni). In this case the variance Vi-j = 1, which 
corresponds to the level of fluctuations in the absence of any correlation between Ni 
and Nj. Variance smaller than one, Vi_j < 1, implies reduction (or squeezing) of 
fluctuations below the shot-noise level and is due to quantum correlation between the 
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Figure 12. Same as in figure [5] except for twice smaller s-wave scattering lengths an 
and ago- 




Figure 13. Illustration of the four regions of the momentum space density, forming 
the quadrants A, B, C, and D on the s-wave scattering sphere, on which we analyse 
the data for relative number squeezing. 

particle number fluctuations in Ni and Nj. Perfect (100%) squeezing of the relative 
number fluctuations corresponds to Vi^j = 0. 

In the context of the present model for the BEC collision experiment and possible 
correlation measurements between atom number fluctuations on diametrically opposite 
sides of the s-wave scattering sphere, we assign the indices i,j = A, B, C, D in Equation 
( l32l) to one of the four quadrants as illustrated in figure [131 The total atom number 
operator Ni in each quadrant T>i within the s-wave scattering sphere is defined after 
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Figure 14. Relative number variance in the diametrically opposite and neighboring 
quadrants, Va^c/b-d and Va-b/c-d, as a function of time. 



elimination of the regions in momentum space occupied by the two colhding condensates 

p r+oo 

N,{t)= / dk^dky dk^h{k,t). (33) 

Jv, J-oo 

Operationally, this is implemented by discarding the data points beyond \kx\ > O.Skr- 
In addition, the quadrants are defined on a 2D plane after integrating the momentum 
distribution along the z-direction, which in turn only takes into account the 3D data 
points satisfying |1 — k'^/k'^\ < 0.28, i.e. lying in the narrow spherical shell k^ ± Sk 
with 6k ^ O.lAkr- The elimination of the inner and outer regions of the halo is done 
to minimize the sampling error in our simulations, since these regions have vanishingly 
small population and produce large noise in the stochastic simulations. 

The choice of the quadrants as above is a particular implementation of the 
procedure of binning, known to result in a stronger correlation signal and larger relative 
number squeezing [Tl| |33] . Due to strong back-to-back pair correlations, we expect the 
relative number fluctuations in the diametrically opposite quadrants to be squeezed, 
Va-c, Vb~d < 1, while the relative number variance in the neighboring quadrants, such 
as Va-b and Vc-d, is expected to be larger than, or equal to, one. The positive- 
P simulations conflrm these expectations and are shown in flgure [HI where we see 
strong (~ 80%) relative number squeezing for the diametrically opposite quadrants, 
Va-c,b-d ~ 1 - 0.8 = 0.2. 

These results assume a uniform detection efficiency off] = 1, whereas if the efficiency 
is less than 100% {rj < 1), then the second term in Equation (l32l) should be multiphed by 
7]. This implies, that for r/ = 0.1 as an example, the above prediction of ~ 80% relative 
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number squeezing will be degraded down to a much smaller but still measurable value of 
~ 8% squeezing {Va-c,b-d — 1 — 0.08 = 0.92). Even with perfect detection efficiency, 
our simulations do not lead to ideal (100%) squeezing. This can be understood in terms 
of a small fraction of collisions that take place with a center-of-mass momentum offset 
that is (nearly) parallel to one of the borders between the quadrants. As a result, the 
respective scattered pairs fail to appear in diametrically opposite quadrants during the 
(finite) propagation time (see also [33j). 

For the symmetric case with (iYj) = (Nj) and (iV/) = (iV?), the variance Vi-j can 
be rewritten as 

V.-, = l + {mgS^-gff], (34) 

(2) 

where the second-order correlation function g^j^ is defined according to 

mm 

Equation flMl) helps to relate the relative number squeezing, Vi^j < 1, to the 

(2) (2) 

violation of the classical Cauchy-Schwartz inequality g\2 > 9ii \ studied extensively 
in quantum optics with photons [3ll[30]. The analysis presented here (see also p3] on 
molecular dissociation) shows that the Cauchy-Schwartz inequality, and its violation, is 
a promising area of study in quantum atom optics as well. 



^If = '-^^7^- (35) 



6. Summary 

An important conclusion that we can draw from the numerical simulations is that 
the predicted widths of the correlation functions are remarkably robust against the 
parameter variations we were able to explore (in Sees. 5.1 through 4). This gives us 
confidence in our physical interpretation of the width as being chiefiy due to the initial 
momentum width of the condensate. The discrepancy with the analytical calculation 
of [25] seems to be primarily due to the different cloud shapes used. The width of the 
halo varies with the parameters we tested in a predictable way and also confirms the 
discussion in Sec. El 

As for comparison with the experiment, the numerically calculated widths of the 
scattering halo and the correlation functions coincide with the experimental ones to 
within better than 20% in most cases. The main discrepancy with the experiment is in 
the ratio of the back to back and coUinear correlation widths. From the experimental 
point of view, these ratios are more significant than the individual widths since some 
sources of uncertainty, such as the number of atoms and the size of the condensates, 
cancel. The discrepancy may mean that the coUinear correlations are not sufficient to 
characterize the size and momentum distribution in the source at this level of accuracy. 
The discrepancies may of course also be due to the numerous experimental imperfections, 
especially the fact that the Raman outcoupling was only 60% efficient, and therefore an 
appreciable trapped rrix = 1 condensate was left behind. This defect may be remedied in 
future experiments. On the other hand, the current simulations neglect the unavoidable 
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interaction of the scattered atoms with unscattered, = condensates as they leave 
the interaction region. This interaction could alter the trajectories of the scattered 
atoms in a minor, but complicated way. Future numerical work must examine this 
possibility further. 

Still, the overall message of this work is that a first principles quantum field theory 
approach can quantitatively account for experimental observations of atomic four-wave 
mixing experiments. This work represents the first time that this sort of numerical 
simulation has been carefully confronted with an experiment. An interesting extension 
would be to examine the regime of stimulated scattering. It has been predicted that 
a highly anisotropic BEC could lead to an anisotropic population of the scattering 
halo [35l [36]- This effect would be a kind of atomic analogue of superradiance observed 
when off-resonant light is shined on a condensate [37l |38]. In addition, our results 
may be useful beyond the cold atom community: theoretical descriptions of correlation 
measurements in heavy ion collisions [39j| may benefit from some of our insights. 
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Appendix A. Duration of the collision 

In order to estimate the collision duration one can consider a simple classical model of 
the collision [22]. Denoting by pi(x, t) and p2(x, t) the density distributions of the two 
condensates, the number of scattered atoms Nsdt) at a given time can be written 



where (Tq = SttAqq is the cross section for a collision of two particles. In this latter 
formula aoo — 5.3 nm is the scattering length between rrix = atoms [T5] . 

The time-dependent density of the two condensates can be calculated from the 
expansion of a condensate in the Thomas- Fermi regime described in |41]. This approach 
suggests two different time scales for the collision duration. First, the separation of 
the two condensates occurs in a time defined by the ratio of the longitudinal size of 
the condensates and their relative velocity tgep = Rx/vr- Taking for Rx the Thomas- 
Fermi radius of the initial condensate, one can show that tsep is on the order of 1 ms. 
At the same time, the condensates expand during their separation on a time scale 
4xp = l/i^s/ = l/f^z — 140 /is. This latter effect appears to be predominant in the 





(A.l) 
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evaluation of Equation (lA.ip and texp can be taken as a definition of the collision duration 
At. The numerical evaluation of Equation (lA.ip gives Nsd^t) ^ 0.66A'"sc(cxd) and the 
estimated total number of scattered atoms corresponds to the experimentally observed 
5% of the initial total number of atoms in the trapped condensate. 



Appendix B. Occupation number of the scattering modes and amplitude of 
the back-to-back correlation 

In order to estimate the occupation number of the scattering modes one needs to 
compare the number of scattered atoms Ngc to the number of scattering modes Nm- 
To achieve this one has to first consider the volume of a scattering mode Vm, given 
by the first-order coherence volume (also dubbed as "phase grain" in [T2| E]). Such 
a volume corresponds in fact to the coherence volume of the source condensate, and 
in practice it can also be deduced from the measurement of the width of the collinear 
correlation function g'^j^{Aki) as one expects in a Hanbury Brown-Twiss experiment. 
For simplicity we match the scattering mode volume Vm to the coherence volume of the 
source condensate in momentum space, 

Vm ^ Pa,{ay,)\ (B.l) 

where (3 is a. geometrical factor which depends on the geometry of the modes. 
Approximating the source condensate in momentum space by a Gaussian oc 
exp[-xV(2a2) - (y2 + z^)/{2al,)], one has (3 = {27rf/^ 

The number of scattering modes Nm can in turn be estimated from the knowledge 
of the total volume of the scattering shell V, 

Nm = (B.2) 

where the volume V is determined from the value of the width of the scattering shell 
St. 

V= j d^k exp[-{k~kry/{25k'^)] (B.3) 

~ A7TV2^kl6k, (B.4) 

for 6k <^ kr- If we apply this estimate to the results of the main numerical example (see 
Sec. 15. ip . we find Nm — 26400. As Ng^ = 1750, this implies an occupation number per 
mode of Nsc/Nm — 0.066. Such an estimate confirms that the system is indeed in the 
spontaneous regime and that bosonic stimulation effects are negligible. 

The simple model of [15] for the back-to-back correlation predicts that its height is 
given by 

gfB{Q) = l+Nm/N,, (B.5) 

Using the above estimate of Nm and the actual value of Ngc found from the numerical 
simulations, we obtain that the height of the back-to-back correlation peak should be 
approximately given by ~ 16. This compares favorably with the actual numerical result 
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of 10.2. Similarly, we obtain the back-to-back correlation peak of: ~ 62 in the example 
with the shorter collision time (compare with the numerical result of 36.6); ~ 18 in the 
example with the smaller collision velocity (compare with 10); and ~ 70 in the example 
with the smaller scattering length (compare with 50). 

Appendix C. Width of the s-wave scattering sphere in the undepleted 
"pump" approximation 

To estimate the width of the halo of scattered atoms beyond the spontaneous regime 
we use the analytic solutions for a uniform system in the so called undepleted "pump" 
approximation in which the number of atoms in the colliding condensates are assumed 
constant. This approximation is applicable to short collision times. Nevertheless, it 
formally describes the regime of stimulated scattering and can be used to estimate the 
width of the s-wave scattering sphere as we show here. 

The problem of BEC collisions in the undepleted "pump" approximation was 
studied in [27]; the solutions for the momentum distribution of the s-wave scattered 
atoms are formally equivalent to those obtained for dissociation of a BEC of molecular 
dimers in the undepleted molecular condensate approximation fiSl 128]. For a uniform 
system with periodic boundary conditions, one has the following analytic solution for 
momentum mode occupation numbers: 



n^it) = =^^sinh^ - AitJ . (C.l) 

Here, the constant g is given by 

g = 2Uopo = , (C.2) 

m 

where Uq = inhaoo/m corresponds to the coupling constant g/hoi [27], and we note that 
the results of [27] contain typographical errors and have to be corrected as follows |42j : 
given the Hamiltonian of (1), with g = Anh'^a/m, the coupling g in (2), (7), (9), and 
(10), as well as in the definition of A(p) after (9), should be replaced by 2g. In the 
problem of molecular dissociation, the constant g corresponds to ^ = Xy/Po [13j, where 
X is the atom-molecule coupling and po is the molecular BEC density. 

The parameter in Equation fIC.ll) corresponds to the energy offset from the 
resonance condition 

HA,^^- (C.3) 
2m 2m ' ^ ^ 

where hkr is the collision momentum; in molecular dissociation, h'^k'^/m corresponds to 

the effective dissociation energy 2/j,|Ae//|, using the notations of [T5] . 

From Equation fIC.ll) we see that modes with — A| > experience Bose 

enhancement and grow exponentially with time, whereas the modes with g'^ — A^ < 

oscillate at the spontaneous noise level. The absolute momenta of the exponentially 

growing modes lie near the resonant momentum hkj., and therefore we can use the 

condition — A| = to define the approximate width of the s-wave scattering sphere. 
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First we write k = kr + A/c and assume for simplicity that kr is large enough so that 

2 

k 



Ak <^ kr- Then the condition g"^ — A^. = can be approximated by 



hkr-Ak^ ^ 



1 - ^ 0. (C.4) 

\ mg J 

This can be solved for Ak and used to define the width Sk = Ak/2 of the s-wave 
scattering sphere as 

6 k mg 47raooPo 



kr 2hkl fc2 



(C.5) 



The reason for defining it as half of Ak is to make 5k closer in definition to the half-width 
at half maximum and to the rms width around kr- 

The above simple analytic estimate (10.50 gives 6k/ kr — 0.05 for the present ^He* 
parameters. For comparison, the actual width of the analytic result f lC.l|) varies between 
6k/ kr ~ 0. 12 and 6k/ kr — 0.027 for durations between 'gt = 1 and gt = 7, corresponding, 
respectively, to t ~ 20 fis and t ~ 140 fis in the present ^He* example. 



Appendix D. Positive-P simulation parameters 

The positive-P simulations in our main numerical example of Sec. [5] are performed on a 
computational lattice with 1400 x 50 x 70 points in the (x, y, 2;)-directions respectively. 
The length of the quantization box along each dimension is L^. = 252 /im, Ly = 20.52 
/im, and Lz = 30.76 fim. The computational lattice in momentum space is reciprocal 
to the position space lattice and has the lattice spacing of Aki = 27r/Lj, giving 
Ak^ = 2.49x10^ m-\ Aky = 3.06xl0^m-\ and AA;^ = 2.04xl0^m-^ The momentum 
cutoffs are fc^""^ = 1.75 x 10^ m~i, fcf^ = 7.66 x 10^ m'^, and k^^""^ = 7.15 x 10^ 
m~^. 

The momentum cutoff in the collision direction, /ci™^^'', is more than 3 times larger 
than the collision momentum kr, and hence it captures all relevant scattering processes 
of interest, including the energy non-conserving scatterings (kr) + (kr) — >■ (3^^) + {—kr) 
and {—kr) + {—kr) {—3kr) + {kr) [H]. In all our figures, the regions of momentum 

are not shown for the clarity of presentation of the main 
halo. These scattering processes, which produce a weak but not negligible signal at 
kx — ±3kr, i.e., outside the main halo are enhanced by Bose stimulation due to the 
large population of the colliding condensate components at k^ — =F^r, respectively. In 
the remaining y and z directions, such processes are absent and therefore the number 
of lattice points and the momentum cutoffs can be smaller. 

Since the momentum distribution of the initial condensate is the narrowest in the 
A;a;-direction, one may question whether the resolution of Ak^ = 2.49 x lO'* m~^ with 1400 
lattice points is sufficient. We check this by repeating the simulations with 4200 x 40 x 40 
lattice points and quantization lengths of = 753 /im and Ly = = 15.4 /im, which 
give smaller lattice spacing Ak^ = 8.24 x 10^ m~^, together with Aky = Akz = 4.08 x 10^ 
m"^, k^^^'' = 1.75 X 10'' m^-*^, and kij^^^^ = k^"^^^^ = 8.16 x 10^ m^-*^. Our results on the 
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new lattice reproduce the previous ones, within the sampUng errors of the stochastic 
simulations. We typically average over 2800 stochastic trajectories, and take 128 time 
steps in the simulations over 25 /xs collision time. A typical simulation of this size takes 
about 100 hours on 7 CPUs running in parallel at 3.6 GHz clock speed. 
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